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Molecular  dynamics  simulation  is  an  important  theoretical  approach  to  many 
problems  in  chemical  physics,  ranging  from  detailed  studies  of  few-body  reaction 
dynamics  and  intramolecular  energy  transfer  [1]  to  the  investigation  equilibrium 
and  nonequilibrium  statistical  mechanics  in  condensed  matter  systems  [2].  A 
simulation  is  essentially  a  numerical  experiment,  however,  and  considerable  data 
reduction  of  the  often  overwhelmingly  detailed  information  generated  by  the 
calculation  must  be  done  before  insight  can  be  extracted.  This  analysis  relies  on  a 
theoretical  viewpoint  to  motivate  the  particular  averaged  quantities  derived  from 
the  dynamical  variables.  For  many-body  condensed  phase  systems  at  or  near 
equilibrium,  statistical  mechanics  provides  the  appropriate  formal  foundation  [3], 
and  the  data  reduction  accompanying  molecular  dynamics  simulations  involves 
calculating  thermodynamic  or  statistical  mechanical  quantities,  such  as  temperature, 
pressure,  density,  pair  correlation  functions,  time  correlation  functions,  and 
transport  coefficients. 

For  few-body  chemical  dynamics  in  the  gas  phase,  recent  advances  in 
nonlinear  science  have  inspired  novel  approaches  to  the  analysis  of  classical 
trajectory  simulations,  and  have  provided  new  insights  into  mechanisms  of  energy 
transfer  and  origins  of  statistical  or  nonstatistical  behavior  [4,5].  Systems  with  a 
small  number  of  modes  are  now  fairly  well  understood,  although  unsolved 
problems  remain  for  systems  with  as  few  as  three  degrees  of  freedom  [5].  Much  less 
is  known  about  dynamical  processes  in  complex  many-body  systems  and  the  nature 
of  nonstatistical  and  mode-specific  energy  transfer  mechanisms  in  large  polyatomic 
molecules,  van  der  Waals  clusters,  and  condensed  phases.  The  details  of  energy 
transfer  in  these  systems  are  being  revealed  modern  ultrafast  optical  experiments  [6], 
stimulating  the  need  for  new  theoretical  approaches  and  insights. 

In  this  brief  report,  we  describe  a  method  for  the  analysis  of  many-body 
molecular  dynamics  simulations,  based  on  a  novel  application  of  joint  time- 
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frequency  approaches  previously  employed  in  nonstationary  signal  processing  [7]. 
The  method  is  based  on  a  generalization  of  correlation  function  methodology  of 
equilibrium  statistical  mechanics  to  the  treatment  of  nonequilibrium  and 
nonstatistical  processes  of  importance  in  ultrafast  phenomena. 

Equilibrium  statistical  mechanics  conceptually  replaces  dynamics  with 
statistics,  and  thus  views  a  complex  many-body  dynamical  problem  as  a  stochastic 
process  with  statistics  characterized  by  time  correlation  functions.  Such  quantities 
can  be  related  to  spectroscopic  and  transport  properties,  and  can  often  be  measured 
directly  by  experiments  [2,3].  An  example  is  the  velocity  autocorrelation  function  of 
a  "tagged"  particle  in  a  many-body  system,  given  by: 

,  U) 

where  the  brackets  indicate  averaging  over  a  statistical  ensemble.  For  ergodic 
systems  at  equilibrium,  the  ensemble  average  can  be  replaced  by  a  time  average 
along  a  single  trajectory: 

T 

c{t)=  \im~  \  Ivis)-  v{t+s)ds 

T  —  1  Jo  ,  (2) 

From  the  perspective  of  signal  processing  [7],  the  dynamics  of  a  many-body  system 
can  be  thought  of  as  a  stationary  random  process;  the  statistics  of  the  process  do  not 
vary  with  time,  thus  allowing  the  time  average  to  be  computed. 

Ultrafast  dynamics  in  condensed  phase  systems  driven  by  chemical  reaction 
or  optical  excitation  differ  fundamentally  from  equilibrium  or  linear  response 
regime  processes.  Initial  conditions  of  extreme  nonequilibrium  nature  are 
generated,  and  the  following  femtosecond  and  picosecond  time  scale  processes  can 
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exhibit  nonstatistical,  nondemocratic,  and  mode-specific  energy  relaxation  pathways. 
This  behavior  cannot  be  thought  of  as  a  stationary  process,  and  indeed,  the 
particular  sequences  of  events  occuring  in  time  and  the  underlying  mechanisms 
governing  them  are  of  central  importance. 

To  generalize  the  idea  of  a  correlation  function  to  ultrafast  dynamical 
processes,  nonstationary  time  series  methods  must  be  employed,  and  concepts  such 
as  evolutionary  spectra  [7,8]  and  time-dependent  frequencies  [5]  must  be  utilized. 
Strictly,  it  is  impossible  to  have  complete  knowledge  of  both  the  time  dependence 
and  frequency  dependence  of  a  signal,  as  the  time-frequency  uncertainty  principle 
guarantees  that  these  two  types  of  information  are  incompatible.  However,  it  is 
possible  to  define  joint  time-frequency  distributions  which  balance  the  lack  of 
certainty  in  a  useful  way.  One  example  is  the  Wigner-Ville  distribution  [7,9].  This 
quantity  is  well-known  in  quantum  mechanics  [9],  where,  with  the  replacement  of 
time  and  frequency  by  position  and  momentum,  a  phase  space-like  distribution, 
known  as  the  Wigner  function,  is  obtained. 

We  define  the  function  0(t,O)),  which  plays  the  role  of  a  time-  and  frequency- 
dependent  effective  "temperature": 


N 

e{t,co)  =  \Y^miei{t,(o) 


(3) 


where  m,-  is  the  mass  of  atom  i,  and  6i(t,Q))  is  given  by. 


0,(f,  (0)  =  \  vSt-  T/2)  •  v,{t  +  T/2)e  d x 


(4) 
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Under  certain  conditions,  and  after  suitable  averaging,  this  quantity  can  be  related 
directly  to  conventional  quantities  in  statistical  mechanics  [3].  By  integrating  over 
frequency,  the  instantaneous  kinetic  energy  is  obtained: 


2x 


J0(f,  £0)d£0=  m.v.(t)  ■v~‘.(t)=  Kit) 


i=l 


(5) 


Alternatively,  for  systems  at  equilibrium,  an  integration  over  time  provides  the 
spectral  density  of  the  mass-averaged  velocity  autocorrelation  function: 


2 

2 


N 

X  m,  I  (v-io)  ■v'^,(t)} 


1=1 


(6) 


Integrating  over  both  frequency  and  time  yields  a  quantity  that  is  proportional  to  the 
time-average  kinetic  energy,  or  temperature,  of  the  system: 


•1  T  '♦’** 

lim-:p-|  0it,  co)dt  do) 


T  -*  -  *  ’0 

N 


i=l 


=  |NfcT 


(7) 


where  k  is  Boltzmann's  constant,  T  is  the  temperature,  and  N  is  the  number  of 
particles. 

The  Wigner  function  &(t,co)  has  the  undesirable  property  of  not  being  positive 
definite,  and  can  experience  rapid  oscillations  as  a  function  of  time  or  frequency. 
Similar  behavior  occurs  for  the  quantum  mechanical  Wigner  function.  A  more 
useful  quantity  is  the  Husimi  distribution  [10],  obtainable  from  the  Wigner  function 
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by  local  gaussian  smoothing.  This  function  is  always  positive,  and  is  simpler  to 
calculate  in  practice.  In  the  example  described  below,  we  shall  employ  the  Husimi 
function  &n(t,ci3),  defined  by: 

,  (8) 

where  cr  is  a  width  parameter  determining  the  relative  smoothing  windows  in 
frequency  and  time. 

To  illustrate  the  method,  we  consider  a  simple  example:  vibrational 
relaxation  of  I2  in  the  van  der  Waals  cluster  I2-Ari3[ll].  Our  model  consists  of  a 
pairwise  additive  potential  function,  with  the  dynamics  of  I2  in  its  B  electronic  state 
represented  by  a  Morse  oscillator  [12]  and  I-Ar  and  Ar-Ar  interactions  given  by  12-6 
Lennard  Jones  potentials  [2].  A  minimum  energy  geometry  of  the  cluster  was 
obtained  by  quenching  a  high  energy  trajectory,  resulting  in  the  I2  impurity 
embedded  in  an  incomplete  solvation  shell  of  argon  atoms.  The  cluster  was  then 
equilibrated  to  an  energy  corresponding  to  a  time-averaged  kinetic  energy  of  3  K.  At 
t  =  0,  the  relative  momentum  of  the  I  atoms  was  boosted  to  correspond  to  the  y  =  20 
quantum  state,  and  Hamilton's  equations  were  integrated  for  300  psec. 

Figure  1  shows  Si{(t,o})  calculated  for  a  trajectory  of  the  I2-Ari3  system.  The 
horizontal  axis  gives  the  frequency  (O  in  wavenumbers,  while  the  vertical  axis  gives 
both  the  value  of  0f{(t,co)  (in  arbitrary  units)  and  the  time  at  which  the  spectrum  is 
calculated  (in  psec.)  Time  increases  with  increasing  vertical  shift  of  the  spectra.  The 
Figure  reveals  a  large  peak  centered  at  100  cm'l  at  the  beginning  of  the  simulation, 
with  very  little  intensity  at  frequencies  below  50  cm'T  The  large  peak  is  due  to  the 
vibrations  of  the  I2  molecule.  It  is  shifted  from  the  harmonic  value  of  128  cm'l  by 
the  anharmonicity  present  at  v  =  20,  which  also  results  in  a  second  small  overtone 
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appearing  at  200  cm"^.  The  initial  absence  of  low  frequency  intensity  is  due  to  the 
lack  of  significant  motion  of  the  cold  argon  atoms.  As  time  progresses,  the  spectrum 
evolves:  (1)  the  frequency  of  the  I2  vibration  increases  with  time,  reaching  about  120 
cm‘l  by  the  end  of  the  simulation,  (2)  the  intensity  of  the  high  I2  frequency 
decreases,  and  (3)  the  intensity  for  low  frequencies  increases  with  time.  These  effects 
arise  due  to  vibrational  relaxation  of  the  initially  excited  I2  molecule,  with  energy 
flow  into  the  low  frequency  interatomic  degrees  of  freedom.  The  I2  frequency 
increases  with  decreasing  diatomic  energy,  due  to  the  anharmonicity  of  the  Morse 
potential. 

In  Fig.  2,  we  investigate  the  dynamics  of  energy  redistribution  in  the  system 
by  defining  two  new  temperature-like  quantities,  corresponding  to  low  and  high 
frequency  components  of  the  motion: 

0  (9) 

.  (10) 

We  take  coc  =  75  cm'^,  a  frequency  that  divides  argon  cluster  and  I2  impurity 
contributions  to  Oh-  These  quantities  measure  the  distribution  of  energy  between 
low  and  high  frequency  degrees  of  freedom.  The  solid  curve  labeled  (a)  shows  the 
time  dependence  of  Khigh,  while  the  long-dashed  curve  (b)  gives  Kiow  In  addition, 
we  plot  the  sum  Ktot  =  f^high  +  f^low  as  the  short-dashed  curve  (c).  Ktot  is 
equivalent  to  a  short  time  average  of  the  instantaneous  kinetic  energy  with  a 
gaussian  weighting  function.  The  total  is  approximately  constant  with  time,  while 
the  decrease  in  Khigh  and  concomitant  increase  in  Kiow  indicates  vibrational  energy 
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relaxation  from  the  initially  excited  I2  into  the  cluster  degrees  of  freedom.  This 
approach  can  be  generalized  to  monitor  the  flow  of  energy  between  "modes"  defined 
in  frequency  space;  this  model-independent  analysis  of  energy  redistribution  has 
advantages  in  highly  nonlinear  systems  such  as  this,  where,  except  for  the  I2 
vibrational  degree  of  freedom,  a  good  set  of  zeroth  order  modes  are  not  known. 

Insight  into  I2  vibrational  energy  relaxation  and  its  effect  on  cluster  structure 
and  dynamics  can  be  gained  from  more  detailed  analysis  of  the  time  and  frequency 
dependence  of  Eq.  (8).  In  Fig.  3  (a)  shows  the  hyperspherical  radius  of  the  system, 
defined  as  p  =  [  +  ...  +  ^5^  ]V2^  where  r,-  is  the  distance  of  the  i**'  particle  from 

the  cluster  center  of  mass.  This  quantity  provides  a  useful  diagnostic,  allowing 
structural  changes  of  the  system  to  be  detected  easily.  Energy  flow  from  the  I2 
impurity  into  the  cluster  induces  isomerizations,  with  the  promotion  of  argon 
atoms  from  their  minimum  potential  sites  across  activation  barriers  to  regions  of 
higher  potential  energy.  Similar  behavior  in  pure  rare  gas  clusters  and  other 
systems  has  been  interpreted  as  a  finite  system  analogue  of  a  solid-liquid  phase 
transition  [13].  This  behavior  appears  as  sudden  increases  in  the  hyperspherical 
radius  of  the  system,  such  as  at  f  =  50, 80,  and  130  psec.  in  Fig.  3  (a).  In  Fig.  3  (b)  we 
show  the  corresponding  time  evolution  of  the  zero  frequency  component  of  &h- 
The  structural  changes  in  the  system  indicated  by  the  time  dependence  of  p( t) 
correlate  well  with  peaks  in  Sif(t,(o=0).  These  peaks  are  caused  by  transient  aperiodic 
motion,  and  occur  when  the  system  is  in  the  process  of  undergoing  a  structural 
change — in  other  words,  when  transition  states  are  being  traversed.  For  early  times, 
the  cluster  modes  are  still  relatively  cold,  and  the  system  migrates  from  one  solid¬ 
like  configuration  to  another.  This  is  reflected  in  the  return  of  GH(t,(o=0)  to  nearly 
zero,  as  seen  for  times  less  that  150  psec.  As  the  argon  degrees  of  freedom  gain 
energy,  &ii(t,(O=0)  experiences  an  overall  increase  with  superimposed  fluctuations. 
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This  correspoiids  to  the  onset  of  diffusive  motion  in  the  now  fluxional,  or  liquid, 
cluster  [13]. 

In  summary,  the  nonstationary  time  series  approach  to  the  analysis  of 
molecular  dynamics  described  here  provides  a  powerful  tool  for  reducing  the  wealth 
of  data  generated  by  computer  simulation  into  an  understandable  and  physically 
intuitive  form.  The  gaussian  smoothed  time- frequency  Wigner  tunction&n  (t,o)), 
which  can  be  interpreted  as  a  time-  and  frequency-dependent  effective  temperature, 
provides  a  clear  view  of  the  dynamical  processes  occuring  vibrational  relaxation  in  a 
model  I2-Ari3  system.  The  approach  provides  a  powerful  tool  for  unraveling  the 
processes  occuring  in  more  complicated  systems,  such  as  large  p>olyatomic 
molecules,  van  der  Waals  clusters,  polymers,  and  solids.  Application  to  these 
problems  will  be  described  elsewhere  [14]. 
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Figure  1 


Figure  2 


Figure  3 


The  time-  and  frequenqr-  dependent  distribution  BH(t>(o)  for  I2 
relaxation  in  Aria-  Horizontal  axis  shows  the  frequency  in 
wavenumbers,  while  the  vertical  axis  gives  both  the  value  of  ©H(t,co) 
and  the  time  at  which  the  spectra  were  calculated. 

High  frequency  (a),  low  frequency  (b),  and  total  integral  of  GuitfOi)  over 
frequency.  See  text  for  discussion. 

(a)  hyperspherical  radius  pit)  vs.  time  for  la-Aria  cluster,  (b)  zero 
frequency  component  of  &H(t,a)) .  See  text  for  details. 
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